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Abstract. In this paper we approximate the radiative transfer equations by 
the method of moments, constructing mesoscopic approximations of arbitrary 
order of the otherwise microscopic system. To define the necessary closure a 
minimum entropy approach is utihzed. While in radiative transfer, the mini- 
mum entropy closure for moment systems up to the first-order moment is well 
known, higher-order minimum entropy closures have not been implemented. 
This is probably due to the fact that the closure cannot be expressed in an- 
alytical form. Our focus thus lies in developing some general results about 
the minimum entropy system and in deriving a numerical closure. By extend- 
ing to higher order, among increasing the precision, we are able to overcome 
difficulties that arise for the first order minimum entropy method. Numeri- 
cal experiments in a 1-dimcnsional domain irradiated by two beams or with 
internal source show the accuracy of this approach. 



1. Introduction 

Radiation therapy has been in use for treatment of cancer for over a hundred 
years and has become one of its most important means. To predict the radiation 
dose, simulations solving an approximation of the radiative transfer equations are 
used. The radiative transfer equations are a system of integro-differential equations 
constituting a balance law for a radiation field in a medium. Currently, most 
applications in clinical use utilize so-called pencil-beam models [T]. However, they 
have fundamental problems handling inhomogeneities of the tissue [ElE], e.g. air 
cavities in lungs or the head, which are extremely vulnerable to irradiation. The 
most accurate methods of calculating the amount of absorbed radiative energy use 
Monte-Carlo models [2]- However, the necessary computation time often exceeds 
the capabilites in clinical environments. 

A third approach is given by the method of moments. Originally developed by 
H. Grad in the context of rarefied gases [9] , through integration over the directional 
variable one averages over all directions, introducing the moments of the particle 
distribution. A transition in scale is achieved: the mesoscopic system becomes 
macroscopic, which is a reasonable simplification, since we are interested in the 
radiation dose and not the exact distribution of all particles. By expanding into 
moments of arbitrary order a hierarchy of moment systems can be defined [6] . A 
closure for the resulting system has to be defined in order to make it solvable. We 
choose minimum entropy closures which enforce that the distribution maximizes 
the physical entropy. The mathematical entropy to be minimized, hence the ti- 
tle, is just the physical one with a minus sign. This approach is inspired by the 
fact that physical systems always tend to increase the entropy. It has become the 
main concept of rational extended thermodynamics [IGj. Jaynes has shown that 
the entropy- minimizing distribution is indeed the most probable one [lOj . Sub- 
sequently, much work in this field has been done by Levermore, who proposed a 
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closure ensuring the hyperbolicity of the system and its abihty to dissipate the en- 
tropy locally [13L I14j . The moment system of order n, closed using the minimum 
entropy principle, is termed Ain- 

Unfortunately, higher order minimum entropy closures are not explicitly express- 
ible, but have to be calculated numerically, which might be the main reason why 
they have not attracted much attention in practice so far. To describe the closure 
the so-called Eddington factor is used. By studying it, we will be able to gain in- 
sight into the system's behaviour, e.g. its ability to handle radiative non-equilibria. 
The one major drawback of the first order minimum entropy method Ali is that it 
cannot handle situations when a net particle flux of zero occurs. As will be shown, 
the M2 model is able to fix this problem. 

The remainder of this paper is organized as follows: in section[2]we introduce the 
underlying transport equations for electrons and photons separately and describe 
the method of moments, after explaining the physical background. In section [3] we 
study the minimum entropy closure in general and the Eddington factor in detail 
to construct the closure of our system. The Eddington factor has to be calculated 
numerically, and we present on e method in section [Jl Section [5] deals with prop- 
erties of different moment systems resulting from the method of moments and a 
minimum entropy closure. Finally, in section [6] numerical results of simulations of 
the transport equations are shown and discussed. 

2. Model Description 

2.1. Physical background. In radiation therapy, beams of protons, electrons or 
photons are used. The common point of interest is radiation dose, i.e. the amount 
of energy deposited in the tissue. Referring to interactions particles undergo within 
the body, we will distinguish between 3 different types: elastic scattering, inelastic 
scattering and absorption. Absorption is the process when a particle loses its energy 
completely. The former type of scattering, elastic scattering, describes processes in 
which a particle's energy remains unchanged, only the direction of flight changes. 
The latter one, inelastic scattering, characterizes the case where the particle loses 
a considerable amount of energy, oftentimes increasing the number of particles 
that contribute to their flux. The probability that a certain particle, interacting 
with an atom, is scattered by a certain angle and gains/loses a certain amount of 
energy is given by the differential scattering cross section, also called scattering 
kernel. As it is in our situation by far more probable that the deflection angle is 
very small, we call the scattering kernel forward-peaked. It is important to check 
whether the model is able to handle this anisotropy. In an opaque medium, the 
radiation is assumed to be near isotropic and macroscopic diffusion models can be 
used [m [IH]. In a transparent medium however, the radiation becomes strongly 
anisotropic where microscopic methods like Monte-Carlo simulations are used. As 
the radiation transfer in our present case belongs to the transitional regime between 
these two extremes, we try to find a model which is macroscopic but nevertheless 
capable of dealing with those non-equilibria. The energy of our particles is the 
relativistic kinetic energy in the case of protons or electrons and E — h ■ v in 
case of photons. A fundamental assumption in this work is that we only consider 
media stationary in space, since otherwise the transfer equation is coupled with 
an equation describing the fluid motion, as particles now would have a prefered 
traveling direction. This is reasonable, since it can be assumed that the patient 
remains stationary while the treatment is taking place. 

2.2. Transport equations. In this section we state the transport equations used 
to describe the transport of electrons or that of photons moving through a medium. 
Let us start by introducing some fundamental notations: let x be the spatial, t 



HIGHER ORDER MINIMUM ENTROPY APPROXIMATIONS IN RADIATIVE TRANSFER 3 



be the time and Q be the directional variable (in a spherical coordinate system): 
= (p, ^/l — fj? cos — ^2 sin (p)^ ■ e describes the energy. It is also possi- 

ble to consider the problem in slab (or planar) geometry: assume the medium has 
infinite length in two dimensions but only finite length in the third one. Imag- 
ine parallel slabs of infinite diameter opposite to each other and a perpendicular 
beam crossing them. That means the (ID-) spatial variable can be identified as the 
penetration depth and the directional variable O collapses to ji, the cosine of the 
angle between direction of flight and the beam. In this case our system becomes 
rotational symmetric. 

2.2.1. Photon transport. We start with the linear Boltzmann transport equation 
(also called radiative transfer equation in this context) 

(1) -dttp{x, n,€,t) + n ■ Vip{x, ft, e, t) = Lptp{x, n, e, t) + Q{x, fJ, e, t), 
where we define the photon scattering operator Lp as 

(2) Lp^{x,n,e,t) := k{x){B{T) - il){x,fl,e,t)) 

+ j C7{x, n' • O, e) • ii){x, n', e, t) dfi' - V(a;, fi, e, t). 

The energy is computed according to e ~ h ■ v with frequency v and the Planck 
constant h. Denote the number of photons by /. Additionally, we define the 
specific intensity ip{x, 51, e, t) = ce f{x, O, e, t), which has the physical interpretation 
that ■ip{x,n,e,t)cos{6)df.dfldAdt describes the amount of radiant energy in the 
interval (e, e -I- de) traveling in time dt through area dA into the element of solid 
angle dfi around fi, where 6 is the angle between 51 and the normal of the area 
dA. The specific intensity describes the radiation field inside the medium. Wc still 
need to define other parameters: the absorption coeflicient, which tells us that a 
photon, traveling a distance ds is absorbed with probability k{x, e) ds. Note that 
absorption and emission share the same coefficient. a{x, 51' • 51, e) is the differential 
scattering cross section describing the probabilty that a photon at point x with 
energy e moving in direction 51' is scattered to direction 51. 

The radiative transfer equation is a balance law for the conservation of energy. 
The left hand side is a transport part which describes how photons travel through 
the medium and the right hand side is a source term accounting for contributions to 
the radiation field. The integral appears due to the scattering of photons into our 
beam (in-scattering) and the last term describes the number of photons scattered 
out of our beam (out-scattering). More specifically, we assume the medium to be in 
local thermodynamic equilibrium (LTE), which allows us to describe the emission 
of photons inside the medium by Planck's distribution function. This is an example 
of a homogeneous and isotropic field, describing the radiation of a perfectly black 
body in thermodynamic equilibrium at temperature T, 

(3) ^ = ^(^)-/£H^)-i)"- 

Here c is the speed of light and k is the Boltzmann constant. 

Keep in mind though, that the radiative transfer equation is in itself only an 
approximate description of the propagation of electro magnetic radiation through 
matter. Two aspects are neglected in this description: firstly, the state of polariza- 
tion of the electro magnetic field is not taken into account and secondly, photons 
are treated as particles, neglecting their wave-like behavior. 
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2.3. Method of moments. In this section we investigate the method of moments, 
which allows us to solve the transfer equation ([T]), by expanding it into a coupled 
system of partial differential equations independent of the angular variable, which 
we will in turn solve numerically. This process can be thought of as averaging over 
all directions. As we will see, this will lead to a closure problem which we will 
tackle in section |31 

A possible motivation for this method is as follows: the transfer equation is a 
mesoscopic equation, between micro- and macroscopic, describing the exact distri- 
bution of all photons in space and time. But in radiation therapy we are in general 
only interested in macroscopic quantities like the energy density. Hence, it makes 
sense to try to achieve this transition in scale by integration. Additionally, we 
can interpret the transfer equation as an infinite system of equations, one for each 
direction, which we want to replace by a finite number of equations. 

From now on, we shall use the following notation 

Definition 2.1. 

oo 

(4) (.) := J J-dQde 

S2 

and introduce the term moment mathematically: we call 

(5) ^^'^ := {fl'^) 
the i*'' moment of ijj, with 

(6) fl°:=l, n^-.^fleind ft' -.^ Q ® . (g) fl. 

Note that is a tensor of i**^ rank and the integration is done component-wise. 
Thus "0^*^ is also a tensor of i'^ rank and an element of r3x -x3^ 

Remark 2.2. The zeroth, first and second angular moment of the distribution are 
called energy density, radiative flux or radiative pressure respectively. 

Remark 2.3. Averaging over all energies leads to so-called grey approximations. 

Definition 2.4. Let us denote by 17^ the linearly independent entries of il*, where 
V + + Now m{Q) is defined as the vector of all fll, Vi < n^k 

(7) m{n) := (1, nl,..., nl nl,..., n-^^.^^f. 

Similarly, we denote by E the vector of the linearly independent entries of the first 
n+1 moments: 

(8) i?:=((V),(f^}^),...,fe_^,„^^^)f. 

In general, to construct a n^^ order moment system {n > 0), one calculates 
all moments of up to order n of every equation, thereby increasing the size of 
the system. That means multiplying the equation by m(f2) and component-wise 
integration, as defined before. The general n*"^ order moment systems for photon 
transport is 

(9) -dt {m{n) ■ tpix, n, e, t)) + V • (r2m(r2) • tpix, n, e, t)) = 

c 

K(x){m{n) ■ {B{T) - i^{x,n,e,t))) 

+{m{n){ J <7{x, n' • 17, e) • ip{x, n\ e, t) dn' - ij{x, n, e, t))) 

S2 
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where the moments of the differential scattering cross section appear in the last 
term, since 

(10) 

1 

{n''{ [ cr(x, n' -ft)- ^{x, n\ e, t) An') = {^t 27r / /^(a;, ^J., e) d^l -ii^ix, n, e, t)). 



=(T('=)(a;,£) 

As can be understood easily, the n}^ order system will always contain the 
(n + 1)''* moment of the distribution, that means first we need to solve the arising 
closure problem, i.e. find . . . The way we choose our 

closure determines the capabilities of our system to model physical situations cor- 
rectly and should additionally guarantee certain desirable mathematical properties, 
e.g. existence of a solution. 



3. Minimum Entropy Closure 

In this section we will introduce the minimum entropy closure and analyze the 
resulting system. The idea is to define the highest order moment to be the respective 
moment of the distribution which minimizes the mathematical entropy of the system 
while reproducing the (given) lower, 0'^ until rt*'\ order moments. We have to bear 
in mind that there is a whole family of distributions which fulfill the latter condition. 
Among those we choose the one which minimizes the entropy of the system, since 
this is the physically most probable one. The mathematical formulation of this 
optimization problem is as follows 

(11) min H{ip) 

s.t. ■ m) — E, 

where H is the entropy, E a vector of prescribed moments, ip and m as before. 
After obtaining the minimizer ipME we set 

(12) ^("+1) {Ci^+'^Pme). 
Definition 3.1. The entropy for bosons is [TUl fT7] 

(13) Hb{^) = {^{{n + 1) ln(n + 1) - nln(ji))) 

where n is the occupation number and it holds that ijj = ^^^-n. 

Additionally, we introduce the following 

Definition 3.2. A^„ denotes the minimum entropy approximation of n"^ degree, 
i.e. whose solution reproduces the first n + 1 moments. 

3.1. Minimum entropy solution. The minizer for equation pT|) is given by [H] 



2hi'^ hv 1 
(14) ipME[a) = -^(exp( — a • to) - 1) , 



where a € R" is the vector of Lagrange multipliers. 

The non-negativity of the solution shows a clear advantage of the minimum 
entropy model compared to the spherical harmonics for example, since in the lat- 
ter unphysical, negative distributions can occur. It is reasonable to simplify the 
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frequency-integral by applying the Stefan-Boltzmann law: 
(15) 

OO 

{^ME{a)) =11 — 5— (exp(— a • m) - l)"MrJdt/ = dstofan / 7 7j df2, 

J J kT J (a - 

S2 52 

where CTstcfan is the Stefan-Boltzmann constant. 

We write our minimizer iPme{<^) in dependency of the Lagrange multipliers, 
instead of the spatial, angular, energy and time variables, as in SI, e, t). The 
Lagrange multipliers have to be determined from the set of constraints {ip ■ m) = E 
and thus depend on the moments i/;'''' < i < n and therefore on (x, i7, e, t) just as 
before. 

Let us state important properties: obviously the underlying distribution is always 
non-negative. Additionally it can be shown, cf. [6], that the moment system closed 
by the minimum entropy closure is symmetrizable hyperbolic. 

Now we take a closer look at the constraints that restrict the distribution -0 and 
its moments or the normalized moments respectively. 

Remark 3.3. From now on, we shall only consider the one-dimensional case, that 
means our model assumes slab-geometry. Effectively f2 is replaced by /i and hence 
all moments become one-dimensional quantities. 

With this assumption it becomes easy to see that the absolute value of i"^ mo- 
ment decreases as i increases, i.e. 



(16) 



Vi e N, 



which tells us that the normalized flux is limited, i.e. 



< 1 which guarantees 



7/j(0) 

that the speed of propagation in our model is limited by the speed of light, a 
fact that does not hold in diffusion models. This follows from a study of the 
characteristic velocities of the system, i.e. the speed of information propagation. 
Now that we calculated the minimum entropy solution, the closure is defined by 

(17) ^("+1) := (M"+VMi=;). 

As mentioned before, the Lagrange multipliers have to be determined from the 
set of constraints {ip ■ m) — E, i.e. it is in general not possible to express 
explicitly in terms of the lower order moments. Therefore we introduce the following 
notation: 

where x is called Eddington factor. One is allowed to work with normalized quan- 
tities since then the Eddington factor no longer explicitly depends on the first 
Lagrange multiplier, see [24j. A further aspect to notice is that the Eddington 
factor is an odd/even function on the odd/even moments, which follows from the 
symmetry of the moments. 

3.2. A4i. It is not possible to derive an explicit formula of the Eddington factor 
for arbitrary dimension n. However, there is an analytical solution for Mi, see [5] 
for a proof: 

V^W, 3 + 4(|^)^ 



V'*"' 5 + 2,/4-3(& 



(19) X( 

V'"' _L . //I _ -^l 

^^(0) , 

See figure [T] for a plot of the (analytically expressible) eigenvalues of the Mi 
system. Let us observe that when the absolute value of the first normalized moment 
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Figure 1: Characteristic velocities of the Mi system 



approaches 1, the distribution has to become more and more peaked around ±1. 
That means the direction of travel for all particles has to be in a smaller and smaller 
cone around the angle arccos(±l). In the limit the distribution becomes a Dirac, 



3.3. A^2- In case of AI2, it is also possible to calculate an (implicit) formula of 
the Eddington factor. The formula is very lengthy and includes certain integral 
expressions and is therefore omitted at this point, see appcndix lA.ll in the appendix. 
To speed up the evaluation we use a rational function with 15 parameters, fitted 
using a least-squares method to the aforementioned formula, see appendix lA.2l and 
[23j . Additionally we will use the fit as a means of verification for our numerical 
computations in the next section. 

By elemental calculus one can show that 



holds. This in turn further confines the domain on which the Eddington factor is 
defined. Corto and Fialkow derived in general realization conditions of trun- 
cated moment problems, but in the present case properties can be derived by direct 
calculations. We want to study what happens when the normalized moments ap- 
proach the boundary of their admissible domain and establish boundary curves, 
which yield the extreme values of the Eddington factor. 

Lemma 3.4. Assume that inequality is sharp, i.e. we are at one boundary of 
the domain of the Eddington factor, the following holds: 

(i) The underlying distribution is a Dirac delta distribution 



(20) 




■0 = S{p — a) 



ae [-1,1]. 



(ii) 



The normalized moments, are of the form 



(21) 



a e [-1, 1] \fi<n+l. 
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(iii) 

^(2) ^(3) ^(1) 

represents a curve and together with equation i21]) defines the admissible 
set of moments for M2- 



Proof, (i) We start by assuming that 
(23) 



holds for some even i > 2. This means that 



(24) / 




but stiU < 1. This equahty can only hold for i/j = ~ a) for some a E [—1, 1], 
as can be seen by using the convolution property of the Dirac delta distribution 

oo 

(25) J f{x)S{x - a)dx = /(a). 

— oo 

(ii) Obvious, using (i) and the convolution property again 

(iii) We can deduce from equation (|20p that when the first normalized moment 
approaches ±1, the second and third normalized moments approach 1 or ±1 resp. 
where the limit for the third normalized moment follows from (j20p . More specifically 

(26) ^ = ±1 ™P^^"' ^^^'^=^^- 

The physical interpretation for either of these two extreme cases is a beam parallel 
or anti parallel to the main axis of interest, i.e. the intensity becomes 5{pL — 1) 
or 5(/i + 1) respectively. Now considering the case of 2 beams and applying the 
superposition principle, we expect the distribution to be a linear combination of 
these two Diracs, i.e. for ^j^^ = 1 we have 

(27) ^ = (^ - i) ■ ^(/^ - 1) + (^ + ^) • + 1), a e [-1, 1] 

and by calculation of moments we get 

from which we deduce the form of the second boundary curve 
(29) ^ = 1 and ^ = ^ = a a^M,!]. 

□ 

The analysis of the boundary helps us to understand the capabilities of the 
model to simulate certain scenarios, e.g. the usage of multiple beams and assists 
us in verifying our numerical computations later on. It will be extended to higher 
orders in section [SJ though the upper bound on the propagation speed remains the 
same. For a more detailed study of different models, but in the electron case, see 
also [5]. 



(28) 
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4. Computation of the Eddington Factor 

As mentioned earlier it is not possible to calculate the Eddington factor for a 
general Mn closure analytically. This section is devoted to the necessary numer- 
ical calculations. Here we present a Monte-Carlo-type ansatz. First, we provide 
Lagrange-multipliers and calculate the Eddington factor as the highest order nor- 
malized moment directly. This is possible since we know the general form of the 
minimizer, 



(30) VmE = T'^Tstefan / ■ . + . . . + )4 

To generate valid moments, observe that the polynomial (q;*^°^ + /ia*^^^ + . . . + 
^"q,("))4^ appearing in the denominator, must not have real roots within the interval 
[— 1, 1], whereas complex roots are allowed. The exception to this is the boundary 
curve. As we know from section 13.31 the distribution tends toward a Dirac, when 
the two normalized moment approach each other, i.e. 

where r is the unique root of the polynomial in the denominator and it holds 
that r e [—1,1]. Technically, in our numerical simulations we circumvent the 
division by zero by calculating the denominator first and if it equals zero, setting 
the integral to 1. The last choice is no approximation but the analytical solution, 
since S{x) da; = 1. It is important to note that the root is only unique for this 
boundary, while for higher orders and other boundaries the distribution becomes a 
linear combination of various Diracs. Hence multiple roots in the interval [—1,1] 
must be allowed. 

Now we randomly, hence titled Monte-Carlo-type approach, choose (valid) 
Lagrange-multipliers a and linearly interpolate the result "^[pio) ^ in the space of 



normalized moments \J^j(oj, ^Jtoji ■ • ■ j ^TotJ- utilize randomness in this ansatz, 
since it is at first glance not clear at all how to transverse the space of normalized 
moments by perturbation of the Lagrange multipliers. For a non-random algo- 
rithm, see |24j . So we simply create enough data points to sufficiently fill our 
domain. Note that figures are shown for 1000 data points for demonstration pur- 
poses, while later computations rely on Eddington factors interpolated from far 
more data points. Interpolation is done using MATLAB, see [50] for details of the 
interpolation algorithm. 

See figure [2] for a comparison between the analytical and the numerical solution 
for Ml, and figure [3] for the M2 closure and its comparison the rational fit. 



5. Properties of the Mn Systems 

Now, with some important properties of the minimum entropy solution, -system 
and Eddington factor under our belt, we shall further investigate the properties 
which can be derived from the hyperbolicity of the system. In order to do so, we 
rewrite the system into the following form: 



(32) 



Ut + F{U)^ = C{U, T) 
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Figure 2: A4i Eddington factor 




with 
(33) 



/^(o)\ 
^(1) * 

^(2) 

V : J 



FiU) 



^(2) 
^(3) 

V : J 



and C{U, T) 



/K(47rB(T)- V(°))\ 
' -(K + a(i))^« 
-(K + a(2))^(2) 



V 



5.1. A^2- First off we want to calculate the characteristic velocities with which the 
system propagates information. They will also allow us to gain insight into how 
the system behaves numerically. The Jacobian matrix J of the system p2p in this 
case, is given by: 



J{U) 



dF{U) 



(34) 



= c 
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Figure 4: Characteristic velocities of the M.2 system 



where we use 



(35) 



(2) 



" I '0(0) 



d 



Since we do not have an analytical expression for the Eddington factor, we use 
symmetric differences to approximate the derivatives. Additionally, for i G {1,2} 



(0) 



dx 



0(0) I 



dx 



d 



\^0(O) ) 



d 



Figure [4] shows the characteristic velocities scaled by the speed of light. As expected 
all eigenvalues are real, since the system is hyperbolic and velocity of the propa- 
gation approaches the speed of light as the normalized moment ^^^^ approaches 1, 
a fact that other models are not capable of simulating. This scenario corresponds 
to the case where all particles are moving in the positive or negative direction of 
the X-axis, depending on the sign of the characteristic velocity. Note especially the 
symmetry of the eigenvalues. 

As mentioned before, it is of great interest to see whether the model can handle 
radiative non-equilibria. Therefore we study how the radiative intensity behaves 
while the angular anisotropy increases. Figured] shows the intensity distribution on 
a grid of normalized moments for different values of /i and the maximum intensity 
for that angle on a logarithmic scale. As would be exptected, the intensity becomes 
more and more peaked for |^| 1. and is isotropic for /i = 0. Figure [S] shows the 
intensity depending on arccos(/i) for fixed normalized moments in polar coordinates. 
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The distribution is highly forward peaked for —> 1 and isotropic for ^^Jj^j = i, 
as stated in section 12.31 while the direction of transportation is perpendicular to 
the X-axis when both normalized moments approach zero which corresponds to the 
second eigenvalue. 




(c) M ^ 1 

Figure 5: Intensity distribution of the A^2 system 



6. Numerical Results 

In this section we discuss numerical simulations using the minimum entropy 
approximation and compare it to other models of radiative transfer. Firstly the 
method of solving those systems will be explained, before we study different sce- 
narios, i.e. test our simulations with different parameters, boundary or initial con- 
ditions. A high order discrete ordinates approximation, cf. [llj . is used as a bench- 
mark in the first test case, while in the second scenario we use an analytical solution 
generated by Su and Olson in [21 . Note that we always assume slab-geometry in 
our simulations, and a 1-dimensional domain. Also, we identify the spatial variable 
X as the penetration depth of the beam. We test minimum entropy approxmations 
of first and second order, denoted by Ali and A^2, against the benchmark. Addi- 
tionally, we add the spherical harmonics methods of order 1 and 3, denoted by Pi 
and P3 respectively. The Pn approach is one of the oldest approximate methods 
for radiative transfer and was first introduced by Eddington in 1926, see [?]• 
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Figure 6: Intensity distribution of the AI2 system; = 0.01 fixed, angles are 
arccos(/^) 
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6.1. Method. As mentioned earlier, the first order system of partial differential 
equations resulting from the minimum entropy method is hyperbolic. We choose 
to use a kinetic scheme |22j . which is a finite volume method, to achieve a spatial 
discretization in the equations. Here, the fluxes over the cell boundaries are com- 
puted, which allows us to prescribe boundary conditions easily. To derive these 
fluxes, incorporating the closure, we use half-moments, i.e. 

1 

(36) := ifi' ■ V)+ := J ^'^ -^dfi {fi' ■ V)- := J f^'-^^dfi 

-1 

where and '-' denote the sign of /i and hence the direction of transport. To 
calculate these half-moments, only minor alterations have to be made to the algo- 
rithm introduced in section 31 Again, computations (and interpolations) are done 
on a grid of normalized moments, as we already have the Lagrange multipliers at 
hand. However, actually only one half-moment has to be computed since we can 
make use of the symmetry in the moments. 

After the discretizations, we end up with a system of ordinary differential equa- 
tions, which in turn is being solved by an adaptive Runge-Kutta method. 

6.2. Simulations. 

6.2.1. First test-case. We set the initial temperature inside the medium to be al- 
most zero and let two rays of light enter at both sides of the domain. The spatial 
domain is a; e [0, 1]. We consider a scattering coefficient of cr = 0.01 and an ab- 
sorption coefflcient of k = 2.5. At the boundary we prescribe the moment of the 
distribution according to tp^^^ = CTstephanT"* with T = lOOOif. 

The test is designed [3] to show the drawback of the otherwise accurate A4i 
model and how A^2 overcomes this problem. 

In the figures we depict the O*** moment of the distribution ip, i.e. the radiative 
energy E, depending on x. See figure [7] for the results. Note that the A4i approx- 
imation produces an unphysical shock and is qualitatively wrong. The A42 model 
fixes this but does not have the same accuracy as the spherical harmonics methods 
Pi and Pg. 

6.2.2. Second test-case. Here the medium is infinite, a: S M. Initially the medium 
is at temperature T — OK for t = 0. A uniform source is applied according to 



(37) Q 



— Xq < X < Xq, < t <to 

otherwise 



with Xq = 0.5 and to = 10. We set the parameters as ct = 0.5 and k = 0.5. 

See figure [5] for a comparison to the benchmark. One can see a distinctive 
dent in all 4 approximations at x = 0.5 where the source is discontinuous. The 
difference between the methods gets smaller as they tend towards zero as the time 
increases. The M2 model considerably improves in accuracy compared to Mi and 
the spherical harmonics approximations. 

7. Conclusion 

Moment methods are an adequate means of treating radiative transfer prob- 
lems. Utilizing those methods it is possible to approximate the transfer (integro- 
differential) equation by a relatively simple system of partial differential equations. 
Still, especially highly anisotropic radiation necessitates a careful choice of the 
closure, as many models have problems when the diffusive regime is left and the 
transition regime is entered. 
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Figure 7: Energy distribution for two entering rays; a = 0.01, k = 2.5, initial 
temperature Tq « OA', boundary temperatures Ti = = lOOOif. 




X 



Figure 8: Energy distribution in an infinite domain a: G M; source in < x < 0.5 
switched on for finite time < t < 10. a = 0.5, k = 0.5, initial temperature 
To « OK. 



Minimum entropy models are able to handle radiative non-equilibria and provide 
good approximations that are computationally cheap. The M2 model corrects the 
Ail model's biggest incapability: its inaptitude to distinguish between radiative 
equilibrium and two identical beams travelling in opposite directions, i.e. it cannot 
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handle a zero net flux. By the employment of the second order method no physi- 
cal shocks occur and the solution always remains positive as opposed to spherical 
harmonics approximations. Also, the propagation of information is always slower 
than the speed of light, a fact that does not hold, for instance, in diffusion models. 

Reasonable next steps would be the extension into more spatial dimensions. 
Then, the moments would cease to be scalar quantities and become tensors. In 
that framework the analysis of the moment methods and the minimum entropy clo- 
sure would be much more involved. Nevertheless, only a realistic geometry would 
allow for modeling inhomogeneities that are non-planar, which are needed to pro- 
vide accurate methods for radiation therapy especially around tissue including air 
cavities such as the lungs and the sinus passage. Furthermore it would be very in- 
teresting to develop new, or modified, algorithms in order to calculate higher order 
minimum entropy closures. This would be supported by a more detailed study of 
the boundary surfaces of higher order Eddington factors. 

Additionally, it would be interesting to study higher order minimum entropy 
methods in different contexts, e. g. plasma physics. 
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Appendix A. Formulae for the M2 Eddington Factor 

A.l. The M2 Eddington factor. 

1 exp(— a) 



63^(erf(i^)-erf(i^)) 
4 b 



exp( 7^^4-^) ((-6)^(8-4a-86) + 2aV-6) 



aexp{a)^/^T{a^ — 66) erf 



la-26\ /I a + 25 



2 J V2 



(38) _2expQ^^^±^±i^) ((-6)i(4 + 2a-46)-2aV-6 
where we defined the error function erf 

X 

(39) erf(a;) := / exp{~t^)dt. 

J 



a and b are the Lagrange-multipliers and have to be determined from the set of 
constraints {ip ■ m) = E, where £^ is a vector of prescribed moments. 

A. 2. Rational fit to the A^2 Eddington factor. A rational fit to equation 
is given by 

, , , ki + k2 - b- k3-b'^ + k4 - b^ + k5 - b^ 
Xfit [a, b)^a ; ; 

(40) +a-^ ■ {l-b)'^" ^ 



3 ^^ ^,^(fc8 + fcg • a^) • (fcio - • 6 + fci2 • 6^ - fci3 • 6^) 



fcl4 + fci5 • 6 

with 
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Parameter 


Value 




0.104107226623813765- 10" 


c 
O 


7 


0.878240820142032813- 10" 


4 


7 


0.834291539331555326- 10" 




7 

K4 


0.182387623676748914- 10" 


5 


7 


0.190298302202491296- 10" 


4 


7 

fce 


0.372238009572563292 - 10" 


4 




0.109376202668751704- 10" 


4 


7 

Kb 


-0.397189030190850628- 10 


3 


7 

Kg 


U.zsy ( oDz ( USi4/ /oizU- iU 


3 


KlO 


-3207.68652563260184 




fell 


2.77365687089333 




A;i2 


2.63485641674103 




fcl3 


0.576015910414080- lO'^ 




A;i4 


-1.17560068853509 




ki5 


0.485606117598845422 





Table 1: Parameters of the rational fit to the A42 Eddington factor 



